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Abstract: Adaptive antenna array processing is widely known to provide significant 
anti-interference capabilities within a Global Navigation Satellite Systems (GNSS) 
receiver. A main challenge in the quest for such receiver architecture has always been the 
computational/processing requirements. Even more demanding would be to try and 
incorporate the flexibility of the Software-Defined Radio (SDR) design philosophy in such 
an implementation. This paper documents a feasible approach to a real-time SDR 
implementation of a beam-steered GNSS receiver and validates its performance. This 
research implements a real-time software receiver on a widely-available x86-based 
multi-core microprocessor to process four-element antenna array data streams sampled 
with 16-bit resolution. The software receiver is capable of 12 channels all-in- view 
Controlled Reception Pattern Antenna (CRPA) array processing capable of rejecting 
multiple interferers. Single Instruction Multiple Data (SIMD) instructions assembly coding 
and multithreaded programming, the key to such an implementation to reduce 
computational complexity, are fully documented within the paper. In conventional antenna 
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array systems, receivers use the geometry of antennas and cable lengths known in advance. 
The documented CRPA implementation is architected to operate without extensive set-up 
and pre-calibration and leverages Space-Time Adaptive Processing (STAP) to provide 
adaptation in both the frequency and space domains. The validation component of the 
paper demonstrates that the developed software receiver operates in real time with live 
Global Positioning System (GPS) and Wide Area Augmentation System (WAAS) LI C/A 
code signal. Further, interference rejection capabilities of the implementation are also 
demonstrated using multiple synthetic interferers which are added to the live data stream. 

Keywords: Global Positioning System (GPS) and Wide Area Augmentation System 
(WAAS) sensor; software-defined radio; controlled reception pattern antenna (CRPA); 
Space-Time Adaptive Processing (STAP); radio frequency interference 



1. Introduction 

Global Navigation Satellite Systems (GNSS) signals are relatively weak and thus vulnerable to 
deliberate or unintentional interference [1]. An electronically-steered antenna array system provides an 
effective approach to mitigate interference by controlling the reception pattern and steering 
beams/nulls. As a result, so called Controlled Reception Pattern Antenna (CRPA) arrays have been 
deployed by organizations such as the U.S. military which seeks high levels of interference rejection. 
However, the increased cost and computational complexity has not been acceptable for the civil 
commercial market thus far. This paper demonstrates a low cost CRPA implementation developed and 
implemented using the flexible software radio approach. 

In the literature, CRPA receivers have been implemented by different approaches. Williams et al. [2] 
used a four-element antenna array and a CRPA system to implement spatial nulling adaptive array. 
The CRPA system used an analogue approach to combine the signals and adopt a correlation feedback 
to derive the weight vector. Konovaltsev et al [3] implemented a MATLAB-based software receiver 
to assess the performance of the beamforming algorithms and found the steering vectors by Direction 
of Arrival (Do A) algorithms. Li et al [4] proposed a beamforming architecture for real Intermediate 
Frequency (IF) signals and calculated optimum weights from given GPS almanac. Heckler et al [5] 
developed a four-element antenna array and front-ends for dual-band L1/L5 signals. They also 
implemented a MATLAB-based software receiver for field test of adaptive beamforming algorithms 
with directions of satellites calculated from DoA estimation. These MATLAB-based software 
receivers do not have real-time capabilities. Arribas et al [6] implemented a real-time hardware and 
software platform which is capable of digital beamforming. They used a Field Programmable Gate 
Array (FPGA) for digital beamforming up to eight antenna elements and delivered the resulting 
spatially-filtered digital signals to a Personal Computer (PC) by Ethernet bus. The signal can be stored for 
post-processing or processed by a real-time software receiver [7,8]. De Lorenzo et al [9] describe various 
approaches for beamforming (pre- and post-correlation) in a post-processing software radio 
architecture. Seo et al [10] use a quad-core Central Processing Unit (CPU) coupled with a Graphics 
Processing Unit (GPU) to implement an all-in-view CRPA software receiver. However, no previous 
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publication has documented the design or implementation of an all-in-view pre-correlation real-time 
CRPA software receiver for GNSS sensors using CPU alone. 

The bit- wise parallel software correlation algorithm [11] has the potential to enable a real-time 
software receiver if digitized GNSS signals have 2-bit resolution. Yet 2-bit resolution does not provide 
sufficient dynamic range when the CRPA design is most valuable — in the presence of high power 
interferer. Thus, the bit-wise parallel software correlation algorithm is not sufficient for any practical 
CRPA implementation which would demand 8 or more bits of sampling resolution. 

This paper completely documents a four-element CRPA software receiver implementation capable 
of processing 16-bit complex samples generating up to twelve beamformed channels for the Global 
Positioning System (GPS) and Wide Area Augmentation System (WAAS) LI Coarse/Acquisition 
(C/A) code signals in real time. 

Conventionally, antenna array system receivers perform CRPA processing with the geometry of the 
antennas and cable lengths known in advance. In this software receiver, the implemented algorithm 
allows CRPA operation without such prior knowledge. Instead, the carrier phase differences between 
satellite tracking channels are used as the constraints of adaptive beamforming algorithm as a 
simplistic model for Minimum Variance Distortionless Response (MVDR) [12] pre-correlation 
beamforming. Further, the implementation includes Space-Time Adaptive Processing (STAP) [13] 
which provides best interference rejection performance but has highest computational complexity 
among interference rejection techniques [14]. 

In order to achieve real-time capability, the most computationally expensive functions are 
programmed in x86 assembly using Single Instruction Multiple Data (SIMD) instructions. We build on 
a structure of SIMD parallel programming [15] and include SIMD assembly code as a real-time 
implementation example since it is the fundamental enabler of such an implementation. Multi-threaded 
programming, again fully documented, is leveraged to fully exploit the multi-core resources of the 
processor [16]. The execution flow is designed to distribute the tasks across multiple cores. For 
reducing the operations of the software correlator, local replicas of code and carrier at zero-phase are 
generated and stored prior to execution. A software technique to correlate the IF data with local replica 
for whole code period without crossing data bit transition is developed and described. 

The implementation is validated using actual IF data streams. The signal processing, positioning 
and beamforming are accomplished by software approach and shown to function in real time. By 
adding synthetic interference to the collected data, the only practical test option given the protected 
nature of the GPS LI frequency band, two test scenarios are used to demonstrate the interference 
rejection performance of the CRPA software receiver. The first interference scenario is with high 
Interference-to- Signal ratio (J/S). The second scenario is with multiple different interferer types from 
different directions. The results are illustrated by angle-frequency responses. 

The key contributions of this paper are summarized as follows: (a) this is the first published work 
which not only demonstrates a real-time all-in- view CRPA software receiver running on a CPU but 
describes the implementation in significant detail with actual code examples critical to the architecture; 
(b) This work leverages a mechanism to perform "blind" beamforming (without prior knowledge or 
calibration). This mechanism is implemented by realizing up to 60 tracking channels (12 tracking 
channels for each antenna element of the four-element array and 12 tracking channels for composite 
beamformed signal), (c) This work provides an approach to estimate Radio Frequency (RF) 
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chain/cabling bias upon startup. This information is not necessary for beamforming in our approach 
but it can be useful to illustrate the functionality of the receiver. Based on this bias information, which 
is used to generate a simplified antenna model, angle-frequency responses are calculated accurately to 
show the functionality of the implemented STAP in the spatial and temporal domains in the presence 
of multiple interferers with different bandwidths. 

The paper proceeds as follows: first, the beamforming algorithm implemented in the software 
receiver is introduced and its mathematical properties are outlined. Then, the mechanism used to build 
the steering vectors without prior calibration is documented. The architecture of the software receiver, 
including hardware and software components, is explained in detail. In order to validate the 
implementation and the interference rejection performance, the software receiver is tested with live 
GPS/WAAS data and synthetic interferers. The results demonstrate the performance of the software 
receiver under multiple interferers in the spatial and temporal domains. 

2. Beamforming Algorithm Used in the Software Receiver 

The primary goal of CRPA is to enhance the carrier-to-noise ratio (C/No) of desired signal and 
reject the interferers. Digital beamforming approaches are used to implement CRPA by combining the 
signals of an antenna array. Various signal combination approaches are well documented in the 
literature, which include frequency-processing, spatial-processing, STAP [13] and Space-Frequency 
Adaptive Processing (SFAP) [17]. Frequency-processing is primarily against narrowband interference. 
Spatial-processing is efficient against both broadband and narrowband interferers, but it can cancel 
only N-l interferers, where TV is the number of elements in the antenna array. STAP or SFAP places 
nulls both in the frequency and spatial domain and can cancel more than N-l interferers if some 
interferers are narrowband. Despite their benefit, the computational complexity of STAP or SFAP is 
high, so implementing them in a real-time software receiver is challenging. Several algorithms have 
been proposed to calculate weights for STAP or SFAP. Some algorithms optimize certain conditions 
with known signal structure of the desired signal such as maximum signal-to-interference ratio, 
minimum mean square error and minimum output power. Other algorithms do not need prior 
knowledge of signal structure but minimize output power to certain constraints such as MVDR [12] 
and Constrained Least Mean-Squares [18]. The constraints can be set to form a beam in the direction 
of satellite or steer a null in the direction of interference. The steering vector associated with the 
direction of satellite can be obtained from either satellite ephemeris and array calibration information 
or even carrier phase differences between elements of an antenna array can be leveraged as a 
simplified model. 

In order to be able to reject multiple interferers with different bandwidths, we adopt the STAP with 
adaptive MVDR beamforming algorithm [19] to implement CRPA shown in Figure 1. After 
down-conversion to zero IF and digitization, the IF signal of each antenna is complex and denoted as 
S n [k]. Each signal is then multiplied by complex weights w n and summed over all antennas to make a 
composite signal as follows: 
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where N is the number of elements in the antenna array. M is the number of time taps. s n [k - m] is the 
signal from n th antenna element and m th time tap. w nm is the weight associated to the n th antenna and m th 
tap. 

Figure 1. Architecture of STAP with adaptive MVDR beamforming algorithm. 
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The MVDR algorithm minimizes the output power and constraints the gain of the direction of desired 
signal to unity as follows [18]: 

minimize \y\ 2 - W H RW 

(2) 

subject to W H C = l 

where R is the covariance matrix of input signals and C is the constraint vector toward the target 
satellite. The covariance matrix R is estimated by computing sample covariance matrix with assuming 
the sample mean is zero as follows: 
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where N s is the number of samples to compute the covariance matrix. The overline sign is defined as 
complex conjugation. For each entity in the R matrix, r ab = r ba a and b stand for the index of antenna 

elements. 

The constraint vector C is composed of a steering vector T as the first to TV* elements and zeros as 
the remaining (M— l)N elements: 



c = 



t N 0 



0 



(M-l)N 



(4) 



T = [t, t 2 ■■■ tj 



Traditionally, the signal of the first antenna element is set as the reference and its component in the 
steering vector is set to 1. The other components are set as phase shifts relative to the reference 
antenna, which are represented by: 



(5) 



where A^is the phase difference based on antenna geometry and the direction of desired satellite. 
Ay n is the phase resulting from the difference of cabling and RF chain including down converter and 
digitization among elements of antenna array. A<p l a can be calculated as: 



(6) 



where p n is the baseline vector of the i th antenna and r l {<p,6) is the unit vector to satellite / shown in 
the Figure 2. 

Figure 2. Antenna geometry and direction of a satellite showing calculation of . 
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The baseline vector p n is known from the array geometry, and r l (<p,6) can be obtained through 
positioning. Only Ay n needs to be re-calibrated whenever any part of hardware of antenna array is 
changed. From the derivation in [18], the solution of Equation (2) is: 



Two matrix inversions are needed to compute in Equation (7). Alternatively, the Frost algorithm [18] 
derived an adaptive approach which iteratively updates the weight: 



where ju is the adaptation step size which can be constant or variable related to covariance matrix. / is 
the iteration count. In our software receiver, ju is calculated by: 



where a is a constant dependent on the number of entries in the R matrix. 
3. Obtaining the Steering Vector without Prior Calibration 

For adaptive MVDR beamforming, obtaining the steering vector is the key to implementing a 
CRPA. However, as mentioned in the previous section, some parameters of the steering vector need to 
be calibrated. Direction vectors toward satellites can be obtained after positioning. A software receiver 
has the flexibility of implementing multiple channels to track multi-antenna signals separately. 
Integrated Carrier Phase (ICP) is one of the tracking outputs of the phase lock loop. The ICP is often 
used to smooth code pseudorange for improving accuracy of positioning. In our software receiver, ICP 
differences between different antenna elements are taken to build the steering vectors instead of 
deriving the vectors from the azimuths/elevations of the satellites and the baseline vectors of antenna 
elements [20]. Figure 3 shows the block diagram of the approach used to obtain the steering vector in 
our software receiver. This approach enables us to perform beamforming without a prior calibration. 

However, if ICP cannot be extracted from individual antenna elements due to weak signal, this 
technique does not work. If the ICP cannot be extracted due to interference, the steering vector is kept 
as previously calculated value. For our implementation, it is assumed that the receiver is initiated in 
interference free conditions. 

The IF data of each antenna is processed by K tracking channels (K = 12 in our case) in which each 
channel is assigned to track one satellite. Each satellite is tracked by N channels (N = 4 in our case) 
assigned to N antenna elements. Once the N channels assigned to a given satellite are simultaneously 
tracked, the ICPs are initialized by filtered phases of the in-phase (I) and quadrature (Q) components of 
correlator outputs. Then, the ICPs are continuously integrated and the phase differences between the 
reference antenna and other antennas are computed at runtime. The differences of ICPs between 
different antennas provide elements of the steering vector. This phase difference, which is the same as 



W = R 1 C(C T R 1 C) 



(7) 




(8) 




(9) 
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Acp l a + Ay iX in Equation (5), can also be used to calculate the cabling/RF chain part Ay n because Acp l a 
is obtained after positioning. 

Figure 3. Mechanism to obtain the steering vectors in the software receiver. 
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4. Architecture of Software Receiver 

The developed CRPA real-time software receiver runs on a PC platform and uses IF datasets as 
input where the detailed description of the data collection hardware is presented in [7]. It currently 
operates on GPS and WAAS LI C/A code signals and supports four-element antenna array. 

There are a total of 60 channels, of which 48 channels are used to process data from four antennas 
(12 channels for each antenna element) and are dedicated to obtain the ICP differences only. The other 
12 channels are used to process the beamformed composite signal using the pre-correlation 
beamforming approach shown in Figure 4. The STAP algorithm, which is capable of rejecting more 
than N-l interferers is implemented in the pre-correlation approach. Each beamformed channel uses 
one set of weights to steer a beam toward a satellite. The weights are updated by adaptive MVDR 
algorithm. The steering vectors, which are inputs to the weight update, are obtained from ICP 
differences between tracking channels as described in the previous section. Hence, the software 
receiver can form beams before positioning in the pre-correlation sense. 

The weight update rate of MVDR algorithm is 1 kHz for adapting highly dynamic interferers. The 
current platform can achieve real-time performance with up to five time taps for STAP processing. 
There exists a tradeoff between the number of antenna elements, tracking channels, time taps, and 
weight update rates. Given the four-element antenna array, our software receiver can reject at most 
three broadband interferers. It can reject more than three interferers if some of them are narrowband 
nature. Positioning is dedicated to the beamformed composite signals. In order to illustrate the 
interference rejection performance, the cabling/RF chain biases need to be calculated after positioning. 
Then, the biases combined with weights are used to calculate the accurate angle-frequency response 
which provides a coarse visualization of the outputs of STAP given the electrical layout of the 
individual elements. 
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Figure 4. Beamforming architecture for steering 12 beams simultaneously using the 
MVDR algorithm. 




4.1. Software Architecture 

The software is developed using Microsoft Visual Studio under 64-bit version of Microsoft 
Windows 7. Most source codes are programmed using C++. The functionalities with high 
computational complexity, such as correlation operation and covariance matrix calculation, are 
programmed by assembly. The component list of the software is in Table 1, and detailed explanations 
of the components are given in the following subsections. 

Table 1. Component list of the software architecture including short descriptions and how 
often each component is called. 



Component 


Description 


Frequency 


Automatic gain control 


Equalize noise power of all elements of 


Call right after reading 




antenna array to full 16-bit representation. 


every 1 ms data from 




Hand-coded assembly. 


solid state drive. 


Software correlator 


Perform 16-bit complex data correlation by 


every 1 ms 




software approach. 






Hand-coded assembly. 




Acquisition/Tracking 


Signal processing and positioning function for 


every 1 ms for 


and Positioning 


GPS/WAAS satellites. 


acquisition/tracking, 






every 100 ms for 






positioning 
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Table 1. Cont. 
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Implement the adaptive MVDR beamforming 
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beamforming 


including covariance matrix calculation and 




-Covariance calculation 


weight updating. 




-Weight updating 


Hand-coded assembly for covariance 






calculation. 




Weight and sum 


Perform the signal combination as Equation 


every 1 ms 
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Hand-coded assembly. 




Bias calculation for 


Perform the calculation of cabling/RF chain 


One time after first 


receiver performance 
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Calculation of angle- 


Perform the calculation of angle-frequency 


every 5 s 


frequency response for 
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performance validation 






System Program 


Schedule threads to achieve real-time 


every 1 ms 




capability. 





Figure 5. Flowchart of the software from reading IF data to validating receiver 
performance. 




4.1.1. Automatic Gain Control 

Software receiver processing starts from reading the IF data from the disk. The MVDR algorithm 
works efficiently when the data received from different antenna gains/RF chain has similar 
distribution. Hence, IF data from each antenna is amplified by a gain term which is managed by an 
automatic gain control function. The objective of the gain control is to equalize the different noise 
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power of all elements of antenna array caused by different RF chains/antenna gains. The gain is 
updated iteratively by finding a maximum value across one ms data as follows: 



MAX = max Ls. 



g = (l Nb -\)lMAX 

G[k + l] = G[k] + j3(g-G[k]) 



(10) 



where s t is / th digitized sample of the IF signal, N s is the number of samples in one ms, Nb is the number 
of bits for representing data, and is updating step. 

4.1.2. Software Correlator 

In the software correlator, the code and carrier of IF signals are wiped off. The 16-bit complex data 
(I/Q) are correlated with local code and carrier replica. Due to its high computational complexity, the 
software correlator is programmed in assembly using SIMD instructions. Based on the SIMD library 
correlator in Hecker and Garrison [15], an assembly code architecture optimized for the specific code 
and carrier tables are created prior to execution. The structures of code, carrier table and input IF data 
are depicted in Figure 6. These tables are formatted as 16-bit short integers. The code table is made by 
N s + 2P S phase for each PRN code where P s is the number of samples between Early-to-Prompt 
spacing and N s is the number of samples in one ms. The carrier table is divided by sine and cosine 
tables and starts from zero phase for each Doppler frequency. I and Q of IF data are bundled together 
to make a pair. The correlation operation is done between IF data with the nearest sample to zero 
carrier phase and local Early/Prompt/Late replicas according to code phase measurement. 

Figure 6. Structure of code table, carrier table, and IF data used in the software correlator. 
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CO: Code 
5: Sine 
CS: Cosine 
/: IF Data I 
Q: IF Data Q 

Ns: # of sample in 1msec 
Ms: # of PRN 
Ps: #of sample in 

E-P spacing 
Fs: # of Frequency 
Format: 16-bit short 
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The correlation operation with complex IF data is performed as follows: 



N.-l 



CJ+jCQ=Z c ° k + A Ih + JQ k ) 

CI = J^CS k x CO k I k -S k x CO k Q h = £C/ t 
CQ = ^CS k x CO k Q k +S k x CO k I k = J^CQ k 



(11) 



k=0 k=0 

where CJ_ and CQ are the real part and imaginary part of correlator output, respectively. The other 
terms are defined in Figure 6. 

Table 2. Procedure to implement software correlator with complex IF data for a single 
sample. 



Step 


Description 


Operations 


1. 


Multiply I and Q of IF data by PRN code and interleave the 


1 multiplications 




result 


2 interleaving 


2. 


Arrange Sine and Cosine tables in complex data format 


4 interleaving 


3. 


Perform complex multiplication by multiply-and-add 


2 multiply-and-add 


4. 


Accumulate the complex results 


3 additions 
1 interleaving 



Using SIMD instructions and 128-bit-wide SSE registers [21], one can simultaneously operate on 
eight 16-bit-wide words by a single instruction. For example, in the step 3 in Table 2, the complex data 
format for executing the parallel multiply-and-add (PMADDWD) instruction is depicted in Figure 7. 
In total, eight multiplications and four additions are performed in one instruction cycle for processing 
two complex IF samples, so it would reduce computational load approximately 12x. 

Figure 7. Complex data format for executing PMADDWD instruction and operations 
inside PMADDWD. 

Bit Count: 128 112 96 80 64 48 32 16 1 



PMADDWD 

(Multiply and Add Packed Integers) 



COiQi 


COxlx 


CO^ 




COoQo 


C0 0 lo 


COoQo 


C0 0 lo 


Ik 


>® 


Ik 




Ik 




Ik 




-Si 






Si 


So 


-CSo 


CSo 


So 


\ 


r \ 


f \ 


f 1 


r 


Cli 


CQi 


Clo 


CQo 



An assembly code of the software correlator, which is compatible to SSE2, SSE3, and SSE4.1 [21], 
is programmed and shown in Figure 8 as an example. This code is built upon the library in Heckler and 
Garrison [15]. Additionally, instructions from the new version of SSE such as interleaving and parallel 
addition are incorporated to further increase performance. This code implements the inner part of the 
summation in Equation (11) and follows the operation procedure in Table 2. Four complex samples are 
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processed in parallel within 30 instructions. This code is the key component to process 60 tracking 
channels in real time and thus its complete description has been included. 

Figure 8. Assembly code example of software correlator using SIMD instructions. 



movupd xmmO, [r8] 

movupd xmml, [r9] 

movupd xmm2 , [rlO] 
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so 
cso 

10 



C02 I COl I COl i coo i coo 



C03 I C03 

Multiply I and Q of IF data by PRN code and interleave result 
C03Q3 | C03I3 | C02Q2 | C02I2 | C01Q1 | COl II | CO0Q0 | CO0I0 



xmm8 — 

xmm5 = C01Q1|C01I1|C01Q1|C01I1|COOQO|COOIO|COOQO|COOIO 

xmm8 = C03Q3|C03I3|C03Q3|C03I3|C02Q2|C02I2|C02Q2|C02I2 

Step 2 . Use Sine and Cosine table to make the complex data format 

xmm6 = S3 | CS3 | S2 | CS2 | SI | CS1 | SO | CSO 

clear xmm7 

move -S to xmml 4 

xmm7 = CS3 | -S3 | CS2 | -S2 | CS1 | -SI | CSO | -SO 
xmm9 = xmm7 



xmm7 
xmm9 
Step 
xmm5 
xmm8 
Step 
xmm5 



SI 
S3 



I 



CS1 
CS3 



I 



SO 
S2 



CSO 
CS2 



CSO 
CS2 



-SO 
-S2 



CS1 | -SI 
CS3 | -S3 

Perform complex multiplication by multiply-and-add 
CI1 | CQ1 | CIO | CQO 

CI3 | CQ3 | CI2 | CQ2 

Accumulate the complex results 
CI3 + CI1 | CQ3 + CQ1 | CI2 + CIO | CQ2 + 
mov xmm5 to xmm6 

move bits 64 . . 127 of xmm5 into 0 . . 63 of xmm6 
xmm5 = CI2 + CIO | CQ2 + CQO 

xmm6 = CI3 + CI1 | CQ3 + CQ1 

accumulate xmm5 to xmml 5 
accumulate xmm6 to xmml 5 



CQO 



In order to avoid performing correlation across data bit transition point, two ms long buffers are 
used in parallel and one ms segment is selected to perform the correlation. The selection approach is to 
use the code phase measurement to set the starting sample of correlation window. Then, this 
zero-code-phase segment is correlated with zero-phase PRN code and carrier tables for whole code 
period. Figure 9(a) shows the architecture of IF data buffer. When new data comes in, the oldest data 
will emerge, the second entry is shifted to the first entry, and the new data is loaded to the second 
entry. Figure 9(b) shows the correlation window which is set by current code phase measurement. 
However, code phase would shift forward or backward due to Doppler frequency. A strategy is 
adopted to account for code phase movement. According to the direction of code phase movement, 
there are three possible cases as shown in Figure 9(c-e) and listed as follows 

1. Forward phase movement to zero (phase Np-lin 1st entry and phase 0 in 2nd entry): perform 
correlation two times for the first and the second buffers. 

2. Backward phase movement to Np-1 (phase 0 in 1st entry and phase Np-1 in 2nd entry): 
perform a correlation for the second buffer and skip next correlation. 
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3. Others: perform a correlation across two entries. 



Figure 9. (a) IF data buffer architecture (b) correlation window (c-e) three cases for 
correlation operation which avoid performing correlation with data bit transition. 



(a) IF Data Buffer Architecture 
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4.1.3. Acquisition/Tracking and Positioning 



These functions adopt open source codes from Greenberg and Ebinuma [22] by replacing its 
interface which controls hardware correlator by our software correlator. Additionally, the WAAS 
signal processing and Viterbi decoder for messages are added to support positioning with the WAAS 
GEO. The structure of carrier tracking loop of software receiver for the specific formats of code and 
carrier table is described in Juang and Chen [23]. And, the phase compensation for carrier wipe-off is 
implemented using a technique in Ledvina et al. [11]. 



4.1.4. ICP Initialization and Calculation of Differential ICP 

The ICP is obtained by integrating IF carrier frequency over time as follows: 

rh k ~ l 
<P P [k] = I codt + (/), - Yj ^['X + 



i=0 



(12) 



where 0 O is the initial phase, co is the Doppler frequency and Ti is the PRN code period which is one 
msfor C/A code. Juang and Chen [23] derived the phase of correlator outputs as: 



0 c [k] = atan2 {CQ[k\, CI[k]) = ^ - 



(13) 



where 0 c [k] is the phase at k ms based on the correlator outputs and the rounded Doppler frequency, cor. 
The correlator outputs, CI[k] and CQ[k], are calculated using the carrier table which has limited 
frequency resolution. co R is the frequency after rounding the Doppler frequency co to the nearest 
frequency grid of the carrier table. Hence, the compensated phase of incoming signal at time tk can be 
obtained by: 
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(p p [k] = 6X k _ l +0 o +6)T I =-0 c [k]+ K R > 1 (14) 

The effect from the rounded frequency in 6 c [k] is compensated in <p p [k] by Equation (14). Since this 
compensated phase is noisy, an averaging filter is used to reduce the noise for estimating initial phase 
A'- 



t^ = ^f ( mod (^ [k-lhaik-llT; ,2k)) + 



( r ^ 



p,[*],2<*<M a (15) 



where 0 p [k] is the averaged phase with less noise, and M a is the number of ms to perform the 
initialization, k denotes k th ms. After the initialization process, the residual of 0 p [M a ] is set as the 
initial phase 0 O as follows: 

fo=modfo p [M a l2x) (16) 

Once the initial phase is obtained, the ICP is calculated every ms by adding a\p^ 1 as in Equation (12). 
In order to build the steering vector, differential ICPs between antennas are calculated for every 
channel. Then, the differential ICPs of several ms are averaged to obtain the steering vectors as follows: 



L 



A^=Zfe[/]-^[/]) 

(17) 

t[ =exp(-yA^ 7 1 ) 

where L is the number of ms to average the differential ICP. / and n stand for / th channel and n th 
antenna. t\ is I th channel and n th antenna element of the steering vector. 

4.1.5. Adaptive MVDR Beamforming 

The adaptive MVDR beamforming function starts from computing the covariance matrix. The 
covariance matrix in Equation (3) is a Hermitian matrix, so only the elements of the upper triangular 
part need to be computed and the remaining elements can be obtained by the conjugation operation. 
The computational complexity of Equation (3) is significant. Indeed, in the four antennas and five time 
taps case, there are 210 elements to be computed. The computational load for each element is N s 
complex 16-bit multiply-and-add operations. With the sampling rate of 4 MHz, 4,000 multiply-and-add 
operations per element need to be completed within 1 ms. Hence, two approaches are adopted to 
expedite the computation. First, multithreaded programming approach is used to assign operations to 
multiple threads for exploiting the resources of multi-core CPU. Second, SIMD assembly 
programming is leveraged to accelerate the complex multiply-and-add operations. The procedure to 
compute covariance matrix initializes with the loading of the complex IF samples from both a th and b th 
antenna element. Then, the multiply-and-add and accumulation operations, as the 3rd and 4th steps in 
Table 2, are performed. After computing the covariance matrix, the adaptation step size is decided by 
Equation (9). Lastly, the adaptive procedure in Equation (8) for each channel is performed to update its 
weight using calculated covariance matrix, steering vector, and adaptation step size. 
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4.1.6. Weight and Sum 

After a set of updated weights is obtained from adaptive MVDR beamforming function, the 
weight- and- sum function performs the combination of IF data as Equation (1). The results are sent to 
software correlator as composite IF data input. Assembly code using SIMD instructions implements 
the weight-and-sum operation by an architecture shown in Figure 10. By arranging the complex data 
format, one multiply-and-add (MADD) instruction processes the complex multiplication of two 
samples and outputs a product including I and Q components. In the end, all products are summed as 
composite IF data. 

Figure 10. Architecture of weight and sum including data formatting and used SIMD 
instructions. 
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4.1.7. Bias Calculation 



Even though this approach can perform beamforming without any calibration/modeling, the RF 
chain/cabling bias information is useful to show the performance of interference rejection of the CRPA 
receiver by angle-frequency response. Hence, these biases are estimated in the software receiver. This 
bias information is not necessary for the beamforming process as previously explained. The bias 
estimation process starts from obtaining the differential ICPs and unit vectors to satellites through 
positioning. The differential ICP for k th satellite between / th and first antenna is given by: 
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&<Pn = 2 ^' r + Aft + N* + e\ (18) 

where A 7 "* is the integer associated to A<p- { , and is the phase error. Then, with known baseline 
vectors, the fractional part of bias can be estimated by: 

/r(A r ,.) = -Lf>od(A^ - l7q, \ r \ ln) (19) 

^ k=l A 

where AT is the number of channels. Only the fractional part of bias is required to calculate the 
angle-frequency response because integer part of bias can be ignored inside exponential function. This 
bias estimation does not consider the biases from mutual coupling and phase-center variations. 

4.1.8. Calculation of Angle-Frequency Response 

After the STAP processing, antenna array response in the direction of (f),0 in Figure 2 and at the 
frequency /is given by: 

N M-l 
n-\ m=0 

where T s is sampling time. There are three dimensions of the response, two for the spatial domain and 
one for the spectral domain. If one would like to show the spatial response such as the gain pattern, 
then the frequency is fixed and gain values of two direction angles are shown on the polar plot. The 
gain pattern is useful to show the spatial performance of beamforming because directional beams and 
nulls are clearly illustrated. As another representation of array response, assuming that one of the 
directional angles of interference is known, one fixes one directional angle and shows the gain of 
the other directional angle versus frequency, which is called angle-frequency response. The 
angle-frequency response is used to simultaneously show the spatial and spectral performance of 
STAP where narrowband notch filter and broadband null steering are simultaneously illustrated. 

4.2. Multi-Threaded Programming and Real-Time Execution Flow 

The system program schedules the threads to execute the software receiver for fully exploiting the 
resource of CPU. Figure 1 1 shows the execution flow plan for the receiver using 4 antennas and 5 time 
taps. To achieve real-time capability, the execution flow must be finished within 1 ms for 1 ms input 
signals. The MultiAnts Thread masters the whole flow to synchronize the threads using event objects. 
At first, MultiAnts Thread sets events to receiver threads and an adaptation thread for initiating the 
individual receiver operation and weight updating. For the composite beamformed signal, RcvrThread 
performs the weight-and-sum using weight from previous iteration, software correlator, 
acquisition/tracking and positioning. The positioning result will be given by this thread. For each 
antenna element, each RcvrWoNavThread performs software correlation and acquisition/tracking only. 
The ICPs are measured by these threads. For adaptation, the AdaptApplebaumThread first sets events 
to multiple CalculateCovarThread to calculate the elements of the covariance matrix, and then 
perform the MVDR algorithm to update the weights. Finally, the MultiAnts Thread uses the differential 
ICPs to build the steering vectors. 



+ Ay n -27tfmT s 



(20) 
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Figure 11. Execution flow of the software receiver within 1 ms. Multithread programming 
is used to make the software capable of real-time execution. 
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5. Interference Rejection Performance of CRPA Software Receiver 



A dataset is collected by the signal collection hardware [24] in open field (low multipath 
environment) to examine the interference rejection performance of the CRPA software receiver. The 
RF front-end of signal collection hardware is a Universal Software Radio Peripheral (USRP) [25] with 
a Bitshark USRP Broadband configurable RF receiver (BURX) daughter board [26]. Its electrical 
details are as follows. Sampling rate is set as 4 MHz. Intermediate frequency bandwidth is set as 
4 MHz. Noise figure is 6.8 dB. Gain is set as 20 dB. Bits per sample at ADC output are 14 bits for both 
I/Q. For injecting the interference, a MATLAB script is written to generate the synthetic interferers 
and combine them with original dataset. The direction and signal specification of interferers can be 
specified in the script. Two scenarios are built. In the first scenario, a single, broadband and high J/S 
interferer is injected. In the second scenario, multiple broadband/narrowband with moderate J/S 
interferers are injected. The power spectral density of broadband and narrowband interferers is shown 
in Figure 12. The 3dB-bandwidth of broadband interferer is 2 MHz and narrowband interferer is a 
continuous sinusoidal wave. 
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Figure 12. Power spectral density of broadband and narrowband interferers. 
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5. 1. Single Interferer with High J/S 

A single interferer specified by 40 dB J/S, 30° elevation angle, and 60° azimuth angle is simulated 
to examine the software receiver performance. The interferer is added into dataset starting from 60th 
second. Figure 13 (left) shows the filtered C/No of all-in- view CRPA. All channels are still in tracking 
when high J/S interferer is injected. Figure 13 (right) shows the C/No of three WAAS GEOs to 
compare performance of CRPA and single antenna. Before injecting the interferer, there is about 5 dB 
gain in C/No of CRPA over single antenna. When the interferer is injected, the C/No of all WAAS 
GEO channels in the single antenna lose lock, but the channels of CRPA continue tracking without 
losing more than 5 dB of the original C/No. Figure 14 shows Earth-North (EN) plots with and without 
interference. The interferer affects the positioning result with worse Circular Error Probable (CEP) 
because the interferer increases the noise level. 



Figure 13. (left) C/N 0 of all-in-view CRPA channels, (right) C/N 0 of three WAAS GEO 
channels (CRPA and single antenna cases). 
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Figure 14. (left) Positioning results of the software receiver without interference, (right) 
Positioning results of the software receiver with interference. 
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5.2. Multiple Interferers 

In order to examine the rejection capabilities of the software receiver to multiple interferers, six 
interferers are injected to the data set as indicated by the sequence shown in Table 3. 



Table 3. Injecting sequence of multiple interferers. 



Sequence 


Interferer Type 


J/S 


IF Frequency 


Direction 


IB 


Broadband 


30 dB 


0.42 MHz 


Azimuth: 0° Elevation: 30° 


2N 


Narrowband 


30 dB 


0.42 MHz 


Azimuth: 60° Elevation: 30° 


3B 


Broadband 


30 dB 


0.42 MHz 


Azimuth: 120° Elevation: 30° 


4N 


Narrowband 


30 dB 


1.42 MHz 


Azimuth: 180° Elevation: 30° 


5N 


Narrowband 


30 dB 


-1.58 MHz 


Azimuth: -120° Elevation: 30° 


6N 


Narrowband 


30 dB 


-0.58 MHz 


Azimuth:-60° Elevation: 30° 



Using a four-element array, the maximum number of broadband interferers that can be nulled using 
spatial adaptation is three. In this scenario, two broadband and four narrowband interferers are injected 
to examine the performance of STAP. Figure 15 shows the sky plots, antenna gain patterns and 
angle-frequency responses of four cases. The antenna gain patterns and angle-frequency responses are 
averaged over all the channels for showing common nulls in the plots. The represented gains are 
obtained by fixing the elevation angle as 30°. 

Each column in Figure 15 stands for each interferer case. The interferers are shown by red dots in 
the sky plot. In the first column, there is no interferer, so no obvious null or notch is present. In the 
second column, when one broadband and one narrowband interferers are injected, one null appear at 
the azimuth angle 0° in the gain pattern and one narrowband notch appear at the frequency of 
0.42 MHz in the angle-frequency response. In the end, total four narrowband and two broadband 
interferers are injected. There are two nulls in the direction of broadband interferers and four notches at 
the frequency of the narrowband interferers. Note that adaptation using the time taps makes notches in 
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the frequency domain instead of steering nulls in the spatial domain. It is the reason why the number of 
interferers which can be rejected is more than the number of elements in the antenna array minus one. 
Figure 16 shows all-in- view CRPA C/No when adding multiple interferers. All channels continue 
tracking throughout the injecting sequence of interferers. This validates the CRPA software receiver is 
capable of mitigating multiple interferers. 

Figure 15. Sky plots and angle-frequency responses of four cases with injected interferers 
marked by red dots in the sky plots. 
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Figure 16. All-in- view CRPA channels C/No under two broadband and four narrowband 
interferers 



All -in -view CRPA C/No 




1B j 2N : 3B : 4N j 5N \ 6N 



0 10 20 30 40 50 60 70 80 90 100 110 120 

Time (s) 



Sensors 2012, 12 



13438 



6. Conclusions 

This paper describes the implementation of a software receiver for GPS/WAAS LI C/A code 
signals to demonstrate the feasibility of CRPA technology for civil applications. The developed 
approach performs pre-correlation beamforming without any prior calibration. The architecture of 
STAP is constructed in a software approach to achieve a CRPA. The optimum weights obtained by 
adaptive MVDR algorithm provide the interference rejection performance capable of rejecting multiple 
interferers. 

The components of software receiver are implemented and constructed efficiently using the SIMD 
assembly and C/C++ multithreaded programming. In particular, the SIMD assembly code has been 
developed and optimized to achieve a fast correlation process of high resolution (16 bits) samples and 
is fully documented within the paper. 

The software receiver implementation is validated, demonstrating it can perform all-in-view pre- 
correlation beamforming up to twelve channels in real-time. The angle-frequency response calculation 
is also included in the software to illustrate the interference rejection performance of the receiver. The 
results of our experiments demonstrate interference rejection capability in the presence of a single 
higher power broadband interference with 40 dB J/S or six lower power broadband/narrowband 
interferers with 30 dB J/S. 

In the next generation of x86 CPUs, a new instruction set called Advanced Vector Extensions 
(AVX) will be incorporated which will have the capacity to further expedite the computation of the 
vector operations. The width of the SIMD register will be 256 bits, two times wider than that of the 
current registers. This would reduce the computational load by two times. Moreover, the AVX 
instruction set introduces a three-operand format where two source operands and one destination 
operand. This would decrease the number of instructions for moving data. Using the AVX instruction 
set, the software receiver would gain the ability to process wider bandwidth signals such as 
GPS/WAAS L5 signals in real time, which require a higher sampling rate. The structure of our 
software receiver is flexible to extend to GPS/WAAS L5 processing with minimal modifications. The 
Federal Aviation Administration (FAA) Alternative Position Navigation and Timing (APNT) study is 
interested in the use of CRPA with the WAAS L5 signal for robust time transfer under GPS 
interference [24]. 
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